An analytical approach to initiation of propagating fronts 
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We consider the problem of initiation of a propagating wave in a one-dimensional excitable fibre. 
In the Zeldovich-Frank-Kamenetsky equation, a.k.a. Nagumo equation, the key role is played by the 
"critical nucleus" solution whose stable manifold is the threshold surface separating initial conditions 
leading to initiation of propagation and to decay. Approximation of this manifold by its tangent 
linear space yields an analytical criterion of initiation which is in a good agreement with direct 
numerical simulations. 



PACS numbers: 87.10.+e, 02.90.+P 

Threshold phenomena are widespread in bistable dis- 
sipative systems. If such a system is spatially extended, 
then fronts switching from one local state to the other 
can propagate. Propagating fronts, or trigger waves, play 
important roles in such diverse physical situations as self- 
heating in metals and superconductors, phase transitions, 
combustion and other chemical reaction waves, and bi- 
ological signalling systems [H [3, H, H, S| to name a 
few. In biology and chemistry they often appear as a fast 
stage of pulse waves in "excitable systems" 0, H, 0, H] ■ 
The question of existence of such waves in particular 
mathematical models is well studied. However, whether 
a propagating wave will actually be observed depends on 
initial conditions. Understanding conditions of initiation 
of propagating fronts or pulses is very important in ap- 
plications. In heart such waves trigger coordinated con- 
traction of the muscle and failure of initiation can cause 
or contribute to serious or fatal medical conditions, or 
render inefficient the work of pacemakers or defibrilla- 
tors 0. In combustion, understanding of initiation is of 
critical importance for safety in storage and transport of 
combustible materials [lo| . 

Mathematically, after the external initiating stimulus 
has finished, the problem is reduced to classification of 
initial conditions that will or will not lead to a propa- 
gating wave solution. This problem is difficult as it is 
essentially non-stationary, spatially extended and non- 
linear, and does not have any helpful symmetries. Yet 
the problem is so important that analytical answers are 
highly desirable even if not very accurate. 

Early attempts of analytical treatment of the initiation 
problems, including the spatially extended ones, used lin- 
ear description supplemented with heuristic conditions to 
represent the threshold [HI, id, lB| and, more re- 
cently, low-dimensional Galerkin styl e ap proximations of 
the partial differential equations 3, 13 

In the last two decades, this problem has been anal- 
ysed from the dynamical systems theory viewpoint [Til . 
HHHSIIlIlli]. These studies identified the impor- 
tance of certain "critical solutions" , whose codimension-1 
(center-)stable manifold acts as the critical surface sep- 
arating the basins of attraction of initiation and decay. 



This understanding was used in sophisticated numerical 
methods of calculating initiation thresholds, e.g. 



21|. 



Here we propose a practical method of defining the 
initiation criteria analytically. The idea is based on the 
linearization of the (center-)stable manifold of the critical 
solution by its linear tangent, the (center-)stable space. 
One would expect that this should work well for initial 
conditions sufficiently close to the critical nucleus. How- 
ever, how close it should be to give a reasonable approxi- 
mation is not clear a priori. We consider a test case with 
very crude initial conditions, in the form of rectangu- 
lar pulses, and the analytical criterion gives surprisingly 
good agreement with direct numerical simulations. 

We consider a one-component reaction-diffusion equa- 
tion 
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with bistable kinetics f{u). As an archetypical example, 
we consider Zeldovich-Frank-Kamenetsky equation sug- 
gested to describe flame propagation [23], which is also 
known as Nagumo equation in its capacity as the fast 
equation in the FitzHugh-Nagumo systern, suggested as 
a simplified model of nerve conduction [2^1, |25j . This 
equation has the kinetics in the form 

f{u)^u{u-e){l-u), 61 = const < 1/2. (2) 

Equations ([1]) have propagating front solutions, 

u — U{x ~ ct A), A = const 

e.g. for 
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We consider a half-infinite cable which is driven away at 
t = from the resting state m = by an instantaneous 
stimulus of amplitude Us and spatial extent Xg at t = 
and/or by a current injection at a; = of amplitude Is 
lasting for time t^, 

ut^u^^ + fiu), ix,t)eR+xR+, (3) 
u.,{0, t) = Isg{t; ts), g(+(X3; t,) = 0, (4) 
u{x,0) — Ush{x;xs), h{+oo;Xs) — 0, (5) 
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FIG. 1: (color online) (a,b) Response to an below- and above-threshold initial perturbation in ZFK equation, (|ll2l8p . Parameter 
values: 6 = 0.13, Is = 0, Xs = 2.10 for both, subthreshold Us — 0.3304831 (a) and superthreshold Us = 0.3304833 (b) cases, 
numerics using central difference centered in space with step hx = 0.15 and forward Euler in time with step ht = 0.01. Dash- 
dotted black lines: initial conditions, bold solid black lines: the critical nuclei, (c) The corresponding critical strength-extent 
curve, separating initiation initial conditions from decay initial conditions, (d) The sketch of a stable manifold of the critical 
solution for the ZFK equation. The critical nucleus is represented by the black dot; the critical trajectories, constituting the 
stable manifold, are shown in black. The family of initial conditions is represented by the dash-dotted line. The bold black line 
is the critical trajectory with initial condition in that family. The sub-threshold trajectories are represented by the blue line, 
while the red lines represent super-threshold trajectories. Note that the point where the initial condition intersect the stable 
manifold is shown as the empty circle. 



or equivalently 

ut = Uxx + f{u) + 2Isg{t-ts)6{x), 

{x,t)eRxR+; (6) 

where S{) is the Dirac deha function (generaUzation for a 
generic stimulus Isg{x,t) is straighforward) . Specifically, 
we consider stimuli of rectangular profiles, 

g{t;t,)^Q{t,-t), h{x;Xs)^Q{xs-x), (8) 

where <d{) is the Heaviside step function. 

Depending on parameters Ms, Xg, Is and tg, problem 
PI4I5P can typically produce either a "decay" solution 
such that maxa; u{x^ t) —> 0^ t oo, see Fig. [IJa), or an 
"initiation" solution such that max^; \u{x, t) — U{x — ct — 
A)| ^ 0, t -> cx) for some A e R, see Fig. ^h). Natu- 
rally, in the even extension (|6l7p , the "initiation" solution 
produces two fronts propagating both ways. Our goal is 
a condition that would predict which of the two outcomes 
will take place for given u^, Xs, Ig and tg. The curve in 
the (tg, Ig) plane, at Ug = 0, separating the two outcomes, 
is widely known as the strength-duration curve. We will 
also consider a similar critical curve in the {xg,Ug) plane 
at Ig = 0, see Fig.[TJc), which we will call strength-extent 
curve. 

We consider first the case Ig = 0, and following [l9| . 
review the fundamental role of the critical nucleus so- 
lution u^{x), which is defined as a nontrivial stationary 
solution of III), i.e. 

+ /("*) ~ 0, ^ const. 



e.g. for Jl]), 

(1 + e)V2 + cosh(a;V2) V2 -59 + 9"^' 

It is then demonstrated that such a solution is unstable. 
Consider linearization of ([1]) near it, u{x,t) — u^{x) + 
v{x,t), v{x,t) <^ 1, then vt = Lu, where h — + 
fu{u^,{x)). Stability of is determined by the spectrum 
of L, 

L0,=A,0,. (9) 

Since L is a Sturni-Liouville operator, all its eigenvalues 
Xj are real. 

Notice that LdxU^ = so u'^(x) is an eigenfunction 
corresponding to eigenvalue 0. By Sturm's oscillation 
theorem, if eigenvalues of the discrete spectrum are or- 
dered so that Ai > A2 > A3 > . . . then egenfunction 
shall have precisely j — 1 zeros. The critical nucleus (x) 
is an even function and has a single maximum at a; = 0, 
so u'^{x) has exactly one zero, and therefore we have 

C2u'^,{x) = (f)2, A2 = 0, 

for some C2 7^ 0. This implies that is unstable, and 
there is exactly one positive eigenvalue, Ai > 0, with a 
corresponding 0i(a;) > 0. The continuous spectrum of L 
is {A} = (-00, Ac], where Ac = lim [9u/(u)] / = 

/'(O) < 0. Hence in the phase space of ©, equilibrium 
is a saddle point, with only one unstable direction. 
Its stable manifold [27] has therefore codimension one 
and, as such, partitions the phase space. One part of 
the phase space corresponds to the decay solutions, and 
the other to the initiation solutions (Fig. Hid)). A one- 
parametric family of initial conditions ([5]), say with a 
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fixed Xs > and the parameter Ug, will cross the stable 
manifold once, say at = For < u*{xs), we 

have decay, and for Ug > u*(xs) initiation. This defines 
the strength-extent curve Us — u*(xs). The role of the 
stable manifold of the critical nucleus the threshold 
surface in the phase space is an empirically verifiable fact: 
it means that the critical nucleus will be observed as a 
transient for any initial conditions sufficiently close to the 
threshold (see Fig.[T](a,b) for t < 100). 

Now we shall use this understanding to construct an 
analytical criterion of initiation. Our idea is to replace 
the stable manifold of by its tangent, i.e. the stable 
space. This implies considering the initiation problem in 
the linear approximation around (x) . Continuing with 
the case Is — 0, we get 

oo 

u{x,t) = + aje^^*(j)j{x), 

i=i 

where for brevity the summation is assumed both 
over the discrete and the continuous spectrum. If we 
choose the eigenfunctions 4>j{x) normalized, then aj = 

J (f>j{x)^{x,0) — it,(x)^ dx. Eigenfunction <j)2{x) — 

u'^{x) is odd, and u{x,0) are even, hence 02 — 0, 

00 

and J2 '^j^^^^4>j{x) ^ as t ^ 00 since Aj < A3 < 

i=3 

for j > 3. Hence in this approximation u{x,t) u^{x) 
if and only if ai = 0. So the equation of the stable 
space, which is an approximation of the critical mani- 
fold, is oi = or 



4>iix) (ush{x; Xs) — u^{x)) dx — 0. (10) 



This is a finite equation for Xs, Us, which provides the 
desired analytical definition of the strength-extent curve. 
For Is ^ 0, we have 



m(x, t) = u^,{x) + ^ Aj{t) cj)j{x), 



J 4)j{x){ u{x^Q) — Uif{x) j dx and 



where Aj{Q) = 

dAj/dt = \jAj -t- 2Is g{t) (l)j{0), which can be solved in 
quadratures for a given g{t), and then the critical condi- 
tion is Ai{+oo) — 0, or 

00 

Ai(0) + 2IsM0) J e-^''9it) dt = 0. (11) 


Now we consider an example with explicit answers. For 
(HI), if 6* < 1, then = 00, and as in [l3|, for u < 6* we 
can approximate 



and then « |0sech^ (^xV9/2^ . In this approximation, 

to solve the eigenvalue problem ([9]), it is convenient to 
change variables v{x) = "tpiz), z = tanh(a:\/^/2), then 

((1-.^)^')'+ (12 -^^1^)^ = 0, ^(±1) = 0, 

solutions of which are associated Legendre functions [2^ . 
In particular, we find that 

Ai = 561/4, 01 (x) Ci sech^(a;v^/2) 

for some Ci 7^ 0. 

For Is — ^ and h{x;Xs) = Q{xs — x), equation (fTU|) 
then gives an explicit equation for the strength-extent 
curve 
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2. , XsVe 

— tanh sech 

TT 2 i I 2 



■ arctan 



n -1 



(13) 



For Us — Q and g{t;ts) = Q{ts — t), we have ^i(O) = 
Ittv^Ci and equation pT|) gives the classical Lapicque- 
Blair-Hill [ll, 12, [l^ equation for the strength-duration 
curve. 



with rheobase 



^rh — 



and chronaxie 



^ MO) ^ 45^^3/2 

Ai /p (l)i{x)u^,{x) dx 64 



(Ai 
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(14) 



(15) 



(16) 



f{u) u{u 



(12) 



Fig. illustrates the quality of the analytical criti- 
cal curves ([T^ and (|14ll5ll6p . both compared to the 
curves obtained by direct numerical simulations for the 
quadratic nonlincarity (fT2|l valid for small 9, and the orig- 
inal cubic nonlincarity ([2|). For the chosen parameter val- 
ues, the error introduced by linear approximation of the 
stable manifold of the critical nucleus is of the same or- 
der of magnitude as the error introduced by the quadratic 
approximation of the nonlincarity. 

In conclusion, we have obtained analytical expressions 
for initiation criteria for a concrete simple example. Such 
criteria were obtained experimentally and numerically 
and any analytical expression was through fitting; we 
have deduced it mathematically ab initio, via a clearly 
defined procedure. The expressions are simple enough 
to be useful in practice, but the procedure of obtaining 
them is probably more important as it can be extended 
to other models. The expression for the strength-extent 
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FIG. 2: (color online) Comparison of analytical predictions 
with numerical simulations, (a) Strength-extent curves for 
rectangular initial conditions, (b) Strength-duration curves 
for point stimulation. Red solid lines: analytical approxima- 
tions, dlSj for (a) and (|14I15I16P for (b). Blue stars ("cub"): 
numerical results for cubic kinetics ([2}. Magenta diamonds 
("quad"): numerical results for quadratic kinetics (|12p. 

curve is specific for the ZFK equation and will have a dif- 
ferent form for a different model. However, the temporal 
strength-duration curve is universal, up to the values of 
two constants, and it coincides precisely with a classi- 
cal form used for over 100 years for analytical fitting of 
empirical data. 

The general principle, linear approximation of the 
(center-)stable manifold of the critical solution, easily 
admits extensions, e.g. for different temporal and spa- 
tial profiles of the initiation stimuli, different initiation 
protocols, possibility of optimization, say with respect to 
the total energy required to initiate a wave etc. 

It also can be extended to other threshold systems, 
whenever the critical solution can be identified, includ- 
ing those having critical solutions which are not critical 
nuclei [111- In such systems, an additional problem is 
anticipated, as one cannot use the even (a: — > —x) em- 
bedding and have to take into account the translational 
symmetry of the problem posed on the whole real axis. 
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